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ABSTRACT 

Galactic nuclei are now generally thought to have density cusps in their centres, and 
the evidence is mounting that as a consequence they are unlikely to be triaxial. Self- 
consistent stellar dynamical models of non-axisymmetric cusps would be an interesting 
counter-argument to this conclusion. We consider 2-d analogues of triaxial cusps: a 
sequence of non-axisymmetric, cuspy discs first described by Sridhar & Touma (1997). 
Scale-free models with potential <i> oc r“ are examined in detail. It is shown analytically 
for 0 < a ^ 0.43 that self-consistent models with positive phase-space density do not 
exist. Numerical solutions of the combined Vlasov and Poisson equations suggest that 
the whole sequence of models with 0 < a < 1 are also unphysical. Along with existing 
work on cusps, we conclude on purely theoretical grounds that galactic nuclei are not 
expected to be triaxial. 

Key words: galaxies: kinematics and dynamics 


1 INTRODUCTION 

The intrinsic shape of elliptical galaxies is a long-standing 
problem, the solution to which may help us to understand 
the process of galaxy formation. Dissipationless collapse and 
merger models predict that ellipticals should generally be 
triaxial with slow figure rotation and anisotropic velocity 
ellipsoids (van Albada 1982, Dubinski & Carlberg 1991, 
Navarro, Frenk & White 1996). But there is no strong obser¬ 
vational evidence for triaxiality (Merritt 1996). Results from 
deprojection of light—fitting stellar kinematics and statis¬ 
tical analysis of the flattening vs. misalignment—are incon¬ 
clusive (Franx, Illingworth & de Zeeuw 1992, Statler 1994). 
In particular, one cannot yet rule out the possibility that all 
ellipticals are either prolate or oblate. 

One purely theoretical approach to the problem is to 
examine whether equilibrium models with strong triaxiality 
exist. Recent observations of the central region of ellipticals 
have found that the luminosity density rises as a power-law 
at small radius (Crane etal. 1993, Gebhardt etal. 1996). 
Schwarzschild (1993) shows that models with strong cusp 
and strong triaxiality cannot be in equilibrium for much 
longer than the diffusion time of the chaotic orbits. This is 
because the regular box orbits are destroyed by the high den¬ 
sity at the centre. The boxlets such as banana and fish orbits 
are not enough to make a self-consistent model without the 
chaotic orbits. Similar study of 2-d non-axisymmetric mod¬ 
els by Kuijken (1993) reaches similar conclusions. Merritt 
& Fridman (1996) discuss 3-d cuspy models which are not 


scale free, and also reach similar conclusions. 

In this paper we examine the existence of cuspy, non- 
axisymmetric models using a new sequence of Stackel models 
by Sridhar & Touma (1997). These models in their simplest 
form are self-similar. They are also 2-dimensional, and as 
such would not make good models for 3-d galactic nuclei. 
Their interest lies in the fact that they are minimally sym¬ 
metric, and hence represent the 2-d analogue of a triaxial 
system. Also very important is the fact that they admit only 
regular orbits, which makes their interpretation much sim¬ 
pler. Sridhar & Touma (1997) postulate that, since banana 
orbits are elongated in the direction supporting the flattened 
potential, self-consistent realizations of the models could ex¬ 
ist. If true, this would be a counter example to the results 
of previous work on non-axisymmetric cusps, and one could 
argue that galactic nuclei may indeed be triaxial. 

After a brief review of the known results of the models 
of Sridhar & Touma (1997) in Section 2, we formulate the 
problem of testing self-consistency in Section 3. Analytical 
results and numerical results are given in Sections 4 and 5, 
and Section 6 concludes. 


2 THE STACKEL POTENTIAL 
2.1 Co-ordinates 

It is convenient to switch between two different coordinate 
systems. In the usual polar coordinates (r, 9), one can define 
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a new parabolic co-ordinate system (^, 77 ), 

^ = r(cos6»-t 1), 
y = r(cos 6 — 1), 

so that curves of constant ^ or constant y are orthogonal 
parabolas. 


2.2 Potential 


In the absence of a central black hole, the self-similar poten¬ 
tial of Sridhar & Touma (1997) is 

ir“[(l-tcos6»)i+“-f (l-cos0)i+“], ae(0,l), (2) 

where in their notation we have set 2K = 1 and 7 = 1 — a. 

This potential is in Stackel form when written in terms 
of the parabolic co-ordinates: 


<E> = 


i+“ _|_ 


(,-n 


(3) 


^ F+{i)-F.{y) 

i-y 

where 

and 

F.{y) = -F+{\y\). 


(4) 


(5), 


( 6 ) 


The potential (2) and the corresponding surface density 
E we write as 


$(r,e)=r“P(0), 

T.{r,e) = r°‘-^s{e). 

where 


(7) 


P{e) = i[(l -f COS0)^+“ + {1- COS0)^+“]. (8) 

The functions S{0) and P{0) are related by Poisson’s equa¬ 
tion, which can be written as 


s{0) = 


00 J7T 

// 


dwd(j) w° 


+ P''{4>) 


[1 -I- — 2w cos{4> — 


(9) 


where primes represent differentiation (see Sridhar & Touma 
1997) and we set 47 r^G = 1. 

As a consequence of (4), the dynamics in the (r, 9) plane 
is analytic. There are two conserved quantities: the energy 
E and a second classical integral 7, given by 

I = 2ipl-iE + F+{0 

= 2ypl-yE F F-{y), 

where {p^^prj) are the conjugate momenta to (^,^ 7 ). The or¬ 
bits of particles are bounded by the co-ordinate lines. The 
shapes of two orbits are shown for illustrative purposes in 
Figure 1. 


9 = 0 



-1.0 -0.5 0.0 0.5 1.0 

r sin 6 

Figure 1. The region covered by two orbits in the model 
with a = 0.5. Both orbits have F = 1; the fat one has 
I = 0.0684 {y/ymax = 0.462), and the thin one has I = 0.1478 
(y/ymax = 0.998). Orbits with negative y are upside-down ver¬ 
sion of these ones: mirror images reflected \n 9 = 7r/2. The zero 
velocity surface E = ^ is shown as a dotted line, on which the 
symbols A and C mark the turning points of the fat orbit; B 
marks the ends of the thin orbit. 


3 TESTING SELF-CONSISTENCY 
3.1 Dynamics 

By Jeans’ Theorem, the distribution function is f(E, I), and 
the associated surface density can be written 

S(r,0) = j /(P,7)dV 


( 11 ) 


The element of phase space in terms of E and I can be 
derived by equating 

dp^dpnd^dy = d^ud^r (12) 

which is true by virtue of the fact that {p^,Pn) are conjugate 
momenta. Hence 

^,^^ld{^y)d{p,p)^EdL 


r d{r,e) d{E,I) 
Now from equation (1) 

d{r,e) 

and from equation ( 10 ) 
d{Pi,Pv) 


2r sin 9 ■ 


1 


diE,I) 2^(7-/_)(/+-7) 

where 

7+ = -Ey + E-{y), 

I-=-E^ + E+{^). 


(13) 


(14) 


(15) 


(16) 


Thus 

where the integration is over the range of {E, 7) for which 
the integrand in real: 7 G [7_, 7+] (7+ > 7_ is implied by 
E > <!>). 
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3.2 Self-similarity 


Self-similarity implies that the distribution function can be 
written in the form 

f = (18) 


Thus, defining 

( 19 ) 

and setting z = ^/E and y = IE~^~^^°‘, the integral relat¬ 
ing E and / (equation 17) can be written 


Q{e) = I z^l--^g{v) 


dzdy 


\/{y - y-){y+ - y) 


( 20 ) 


where y± are defined in the Appendix. 

The region of integration in equation (20) is that 
bounded by the curves y = y±{z, 9). It is illustrated in Fig¬ 
ure 2. 

The maximum value of \y\ is given by the turning point 
of y+ with respect to variations in « and 9. It is easy to show 
by differentiation that the maximum is 
a 

= (1 + a)i-t-i/^ ■ (21) 


The orbits of particles with y = ±ymax are along the 
co-ordinate lines, and hence are thin bananas (cf Bishop 
1987)—parabolic thin shells with their tips lying on the sur¬ 
face i? = Note that the largest value of y in the shaded 
region in Figure 2 is less than j/max, and hence there is no 
thin banana with y > 0 for a — 0.5, 9 — tt/S. The smallest 
value of y is however —2/max, so there is a thin banana with 
y <0. 


3.3 Self-consistency 

The problem we would like to solve is to find a distribution 
function which self-consistently generates the density and 
the potential in equation (7); i.e. we seek g{y) such that S{9) 
is given by equations (20) and (19), with Poisson’s equation 
(9) simultaneously satisfied. 

Since the model density is symmetric about 9 = 0 and 
7 r/2, we should always weight equally a banana orbit and its 
reflection in 6^ = tv/2; thus g{y) is even, and we need only 
consider positive y. We define K+ as the contribution to Q 
from a single positive value of y. 

K+{9,y)= f ^1/^-2 -^^ ==. (22) 

J V (y - y-)iy+- y) 

The corresponding quantity K- for negative y is just the 
reflection of K+ in the major axis, 9 = tv/2: 

K-i9,y) = K+{n-9,-y). (23) 

For each value of 9 there is a contribution to Q from ±y, 
and we write 

Q{9)= J K{9,y)g{y)dy, (24) 

y>0 

where 

K(9,y) = K+{9,y) + K.{9,-y). (25) 



Figure 2. The region of integration in equation (20) for a = 0.5, 
6 = 7r/3 is bounded by the curved lines. For g{y) piecewise con¬ 
stant it is convenient to split the integral into regions as shown 
by the boxes (here yj = 0.06 and Sj = 0.02). 



Figure 3. The function K+{0,y) for the two values of y corre¬ 
sponding to the orbits in Figure 1. The range of 9 has been ex¬ 
tended to [0, vr] for clarity. Thus, by equation (23), K— can be read 
off from the 9 > 7r/2. The solid curve with the most pronounced 
spike is the thin orbit with j//j/max = 0.998; the dotted curve 
with a lower, broader spike is the fat orbit with j//j/max = 0.462. 
The symbols A and C mark the turning points of the fat orbit; 
B marks the ends of the thin orbit (cf. Figure 1). 


The function K{9,y) is illustrated in Figure 3 for two 
values of y. The banana shapes corresponding to these orbits 
are shown in Figure 1. The discontinuity and the cusp {A 
and C) correspond to the the top corner (A) and the bottom 
corner (C) of the fat banana in Figure 1. The end of the 
thin banana gives rise to similar features {B, cf. Figure 1) 
which are closer together. Since neither orbit is very sharply 
bent (cf. Figure 1), there is no contribution to K at small 
9. In fact the only orbits which come near 9 = 0 (the minor 
axis) are fat ones—those with small y. The thinner orbit in 
Figure 1 does not contribute to the density at small 9. This 
illustrates that the thin shells cannot be arbitrarily sharply 
bent; i.e. all thin bananas subtend a finite angle. 
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3.4 Matrix formulation 

Unfortunately, y± depend on 2 in such a non-trivial way that 
the integral over 2 in equation (20) cannot be simplihed 
further analytically. But the integral over y can be done 
analytically for constant g, so our approximation will be to 
make a grid in y and make g(y) piecewise constant on the 
grid; thus we set g = gj = constant for y € {yj ± Sj/2). We 
also need a set of values of 6 = 6i at which to try and match 
the value of Q with that produced by our guess at g. Thus 
we recast equation (24) as a matrix equation: 

Qi ~ 'y ' Kjj9 3 (26) 

3 


with 


yj+Sj/2 

Kij= J K{ei,y)dy. 

yj 


(27) 


We obtain 


Kij = dz z 


ljct-2 


%j 

{sin~^ h{yj + (5j/2) - sin"^ h{yj - 5^/2) 

+ sin~^ h(—yj + (5j/2) — sin”^ h{—yj — 5j/2)} , 


(29) 


where 


1 , 

y < y+\ 

h{y) = 

2 (j/-j/-) 

- 1, 

y- <y < y+; 

1 

1 

+ 

-1, 

1 

lA 



with y± evaluated at {z,0i). With h defined this way, the 
limits of integration in equation (28) can be taken to be 
(0,1), but it is more efficient numerically to split the interval 
into three parts according to the intersections of the strip 
y e {yj ± Sj/2) with y±. An example is shown in Figure 2: 
the integration is carried out over each of the three boxes 
and the results added together. 

To compute Q{6) (equation 19), we need to compute the 
angular part of the density, S{6). For our purposes it will 
suffice to expand S and P in a Fourier series and truncate 
at high order: 


n 

sio)^ E Sk cos{k6) (30) 

fc =0 

n 

p{e) = E PkCos{k6) (31) 

fc =0 

where Pk can be evaluated trivially from equations (2) and 
(7) by standard integrals, and 

= (32) 

(Kalnajs 1971, Syer and Tremaine 1996). Since P{0) is not 
only even, but also symmetric about 0 = n/2, pk and Sk 
vanish when k is odd. We truncate at fc = 30, which leads 
to negligible error for a > 0.5. 


4 ANALYTIC RESULTS 

A necessary condition for self-consistency is that the orbits 
can produce the flattening of the model. An indicator of the 
flattening is the ratio of Q on the minor axis to that on the 
major axis at the same radius. Namely 


Q(o) 

- Q(V2)- 


(33) 


We can also define a similar quantity for each orbit y. 


^ K{f),y) 

A(7r/2,y)- 


(34) 


To make a self-consistent model, the value of C for the den¬ 
sity must be within the range of the values given by the 
orbits: i.e. 


min[Cji] < C{a) < max[Cji], 


(35) 


for a self-consistent model with a potential cusp a. {First 
mean value theorem, inequality 12.111, Gradshteyn & Ry- 
zhik 1980.) 

The maximum of Cy is set by the fattest banana orbit 
(y —> 0), because it spends most of its time near apocentre, 
near the minor axis (cf. Figure 1). For such orbits, 

{ 00, a<l, 

^ 1 (36) 


The minimum of the Cy occurs for the thinnest (closed) 
banana orbit {y i/max). For such orbit, Cy = min[Cy] = 
l/\/2. (See the Appendix for derivations of these limiting 
values of Cy.) 

Figure 4 shows C{a) as a function of a. The inequality 


V 2 


< C{a) < max/Cy] 


(37) 


holds only for a ^ 0.43. So models with 0 < a 0.43 cannot 
be made self-consistent. The inequality (35) is only a nec¬ 
essary condition for self-consistency, so models with higher 
a may or may not be self-consistent. Note, however, that 
max[C'y] ^ 1 and C ^ 1 as a ^ 1, so (35) is only marginally 
satisfied. This may be a hint that it will be difficult to find 
self-consistent models for a ~ 1. To draw stronger conclu¬ 
sions we must try to solve equation (26) numerically. 

To conclude the section on analytic results, we aban¬ 
don for a moment the goal of self-consistency, and remark 
that when g{y) = constant, Q is independent of 6. Thus 
E oc and the contours of surface density coincide 

with the contours of the potential. As a approaches unity, 
E approaches a constant while 'F remains highly flattened. 
Intriguingly when a = 1, g{y) — constant is the trivial self- 
consistent solution. This model is however unphysical: E is 
constant but it vanishes for finite <F. 


5 NUMERICAL SOLUTIONS 

The problem of finding self-consistent solution reduces to 
solving equation (26). A minimum requirement for a physi¬ 
cal solution is that the solution g{y) be positive definite. It 
should also be reasonably smooth. 
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Figure 4. The ratio of the density at the minor axis vs that 
at the major axis for the Sridhar Sz Touma (1997) models. The 
shaded regions indicate those forbidden from self-consistency by 
inequality (35). Models with C < l/y/2 (0 < a ^ 0.43) cannot be 
made self-consistent. 


10-2 


10-3 




10 '^ 


10 - 6 ^1 , I , I , I , I , I , I 

0.4 0.6 0.8 1.0 

a. 

Figure 5. The residuals in the form of of best fit orbit models 
as a function of a. Triangles are for Nq = 11, and circles for 
Nq = 20. From top to bottom in each case Ny = [25, 50,100]. All 
models have significant residual (y^ ^ 10“^^). 


Thus we minimize the function 

i=l 


and require 


The smoothness is constrained through 

S[g] = (gi - g2f + {gNy - gNy-if 

Ny-l 

+ T! ~gj-i -gj+if 

3=2 


(38) 


(39) 


(40) 



Figure 6. The target density Q (solid line) in the case 
a = 0.9. The dashed line shows the best fit model with 
(NQ^Ny) = (50,100). The dotted line shows the best fit model 
with {NQ,Ny) = (50,200). The corresponding values of s,re 
0.0026 and 0.0023 respectively. In both cases A = 10“^^, and 
consequently g(y) is very far from smooth. 

and a small tunable parameter A. We set the characteristic 
scale for the deviation of the solution from the target to be 

'^ = Q(f)-Q(o). (41) 

This problem is solvable with standard routines such as 
QPROG of IMSL and E04NAF of NAG. A source code is in 
Hanson & Lawson (1995). Also an excellent general discus¬ 
sion of the method is in Press et al. (1992). We use a uniform 
grid in 9 with Ne = [11,20,50], and a logarithmic grid in 
y with Ny = [25, 50,100, 200]. We apply a minimal smooth¬ 
ing constraint, A = 10“^^ which leads to solutions with a 
large amount of structure on the smallest scales. Taken lit¬ 
erally, this kind of distribution function is almost certainly 
unphysical, but it gives the most generous constraints on 
self-consistency; if we cannot find even highly oscillatory so¬ 
lutions, then smooth ones probably do not exist. We would 
consider models to be self-consistent if was consistent 
with zero plus rounding error: i.e. ^ 10“^^. 

Figure 5 shows the value of for a range of models. 
Figure 4 suggests that self-consistency should be least con¬ 
strained where the solid line is furthest from the forbidden 
regions. Hence we see in Figure 5 that the smallest value of 
occurs for a ~ 0.7. Figure 5 also shows that additional 
resolution at a = 0.7 gives larger values of and thus pro¬ 
vides a tighter constraint on self-consistency. None of the 
models has x^ consistent with zero. 

We note that a more stringent test of self-consistency 
would perhaps be the maximum residual Amax (as opposed 
to the mean square residual y^). From Figure 6 we see that 
typically the residuals cluster near one or both of the axes 
{6 — 0, tt/2). Thus Amax can be much greater than ■ 

5.1 Fricke expansion 

To gain some insight to the source of non-self-consistency, 
we introduce a basis function approach to solving the Vlasov 
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Figure 7. The target density Q (solid line) for a = 0.99. Dotted 
lines are the contribution to Q from a distribution function with 
g(y) = 3 /", for n = [0, 2, 4, 6]. 


equation. Consider the distribution function obtained by set¬ 
ting g{y) = y" with n an integer. This is the equivalent of 
the method of expansion of a 3-d axisymmetric system due 
to Fricke (1951). Thus we set g{y) = y" in equation (20) 
and define 


Q„(e) 


l/a-2 n dzdy 

z ' y — 

y [y - y-)(y+ - y) 


(42) 


If g{y) is developed in a Taylor expansion about y = 0: 


g{y) = a^y’^. 

n 


(43) 


then Q{6) is a Fricke expansion: 

Q{e) = a„Qn{e). (44) 

n 


An expression for Qn{0) is derived in the Appendix. Q„{9) 
is plotted for n = [0,2,4, 6], a = 0.99 in Figure 7. All the 
Qn with n > 0 have broadly similar shapes. They are flat 
for 6^0 and have negative curvature; they rise and reach 
a maximum at around 6 = O.drr; and they are all flat again 
a,t 9 — 7r/2. The maximum comes from adding up a series of 
cusps like the ones in Figure 3. Its position is fixed roughly 
by the minimum angle subtended by the thin bananas. The 
shape of the function Q, although it is very close to being 
constant, is rather different. We can understand broadly why 
the residuals in Figure 6 cluster at the axes. Given this, it 
looks rather unlikely that a Fricke expansion (44) of Q would 
converge. 

To find the best fit Fricke expansion for a given model 
one can minimize the deviation of equation (44) from the 
true value of Q{9) with respect to variations in a„. In gen¬ 
eral it is difficult to control the values of a„ so as to keep 
g{y) positive definite. For this reason we do not recommend 
Fricke expansion as a method of solving the Vlasov equation. 

For a ~ 1 we know that g(y) is approximately constant, 
so there is some hope that the a„ will be small for n > 0, and 
hence that g{y) will remain positive. We tried fitting Q{9) 


for a = 0.99 with a Fricke expansion. Very poor fits result up 
to n = 14, and are improved very little by going to n = 22. 
The corresponding g{y) turns out not to be positive definite, 
and it is highly oscillatory (more so for the higher order fit, 
showing the non-convergent nature of the expansion). 


6 CONCLUSIONS 

We attempt to build self-consistent models for the sequence 
of Stackel potentials of Sridhar & Touma (1997). Our nu¬ 
merical solutions suggest that such models do not exist for 
the whole range of 0 < a < 1. We also found a simple an¬ 
alytic argument to show that models with 0 < a 0.43 
cannot be self-consistent: all orbits put too much mass near 
the minor axis of the potential, and too little mass near the 
major axis. Unfortunately the argument does not apply to 
higher a. 

We only considered the scale-free models of Sridhar 
& Touma (1997). There are two important generalizations 
which are not scale free: those with central black holes, 
and those which are linear superpositions of models with 
a e (0,1). Consider first a black hole model consisting of a 
model with a € (0,1) but including a black hole of mass M. 
(The black hole component can be thought of as a model 
with a = —1.) Close to the black hole there are lens or¬ 
bits in addition to the bananas, and there is a chance that 
a model composed largely of lenses could be self-consistent 
(Sridhar & Touma 1997). However, at large radius the black 
hole is insignificant, and the system should approach the 
scale-free case. To be more precise, as r ^ oo the only or¬ 
bits which pass close to the black hole will be increasingly 
fat bananas (y —> 0). In the scale-free case these orbits have 
a density biased heavily toward the turning points at the mi¬ 
nor axis, and are therefore unlikely to be useful for building 
self-consistent models even in the presence of a black hole. 
Thus we expect self-consistent models to face the same dif¬ 
ficulties as in the present work. Simple linear superpositions 
of a finite number of models with a € (0, 1) will have similar 
asymptotic problems. 

Our results show that regular orbits in a particular set 
of cuspy non-axisymmetric potentials do not have enough 
freedom to make a self-consistent model. The models we 
study are 2-dimensional, but because of their low symmetry 
they are analogous to triaxial systems. Had self-consistent 
models been found it would have been an interesting counter 
argument against the conclusions of existing work on triaxial 
cusps and galactic nuclei. In fact, we are in agreeement with 
the conclusions of Merritt & Fridman (1996) with respect to 
a rather different potential set. In this case the conclusion 
was that the box orbits required to build non-axisymmetric 
models are made stochastic by the central cusp. In contrast, 
the potentials we use do not admit stochastic orbits, and yet 
still no self-consistent equilibria were found. The existence 
of triaxial galactic nuclei is again thrown into doubt by this 
result. 
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APPENDIX 
Definition of y± 

The boundary of the region of integration in equation (20) 
is defined by the quantities 

V- = —uz ' [1 — u z), 

u (45) 

y+ = vz - v°‘z), 
where 

_ 2 _ 

{u,v) = {e\v\)^ “ 

1 

= (1 +COs6f, 1 — COs6f) P{9)^a. 

From equations (2) & (7) it follows that P{9) > P{n/2) = 1, 
u{9) > u{0) = u{tv/2) = 1 and v £ (0,1). 


Limiting values of Cy 

Next we derive the value of Cy (equation 34) for the limiting 
cases J/ ^ 0 and y ymax- 

At 0 = 0 we have (m, w) = (1,0) and = 0 while 
—y- = x{z) where 

x{z) = z^^°'{l — z). (47) 

Thus K+ = 0 and 
K(0,y) = K.{0, -y) 

_ 1 f ..1/0-2 dz ( 48 ) 

Vv J V^-y' 


At 9 = -K12 we have u = v = 1 and y+ = 
Thus K- — K+ and 


K{'xl2,y) = 2K+{n/2,y) 


= 2 


^l/a-2 


dz 


-y- = X. 


(49) 


y/{x - y){x + y) 


Asy ^ 0 the numerator in (49) diverges at z = y. Thus 
the value of the integral is dominated by the contribution 
from 2 ~ 0, and we can replace x by 2 ^^“. Changing variable 
to ui = z^^°'ly, we obtain 

K{-K/2,y) ^ 2ay~°‘ J (50) 

whence, extending the limits of integration to [1, oo), 

K{'x/2,y) ^ y~°‘aV^ _^ as y ^ 0. (51) 

This expression is valid for a > 0, but for a ~ 0 it is not very 
accurate. This can be rectified by adding the next highest 
order (logarithmic) term: 


K{-K/2,y) ^ ay +2{1 +a) In I (52) 

1(2) y 

For a < y as J/ ^ 0, the integral in equation (48) 
converges, so: 


K{0,y) 




as y ^ 0. 


(53) 


For a > y the same trick of replacing x by z^^°‘ works, and 
we obtain 


Ki0,y) 


— a 

y a 


F(a) 


Thus 


Cy 


y _^ V 2 a 1 V 2 ^ 2 

rwr(f) ’ 


a < i, 
a > i, 


(54) 


(55) 


as y ^ 0. This is the maximum value of Cy. 

As y —> ymax, x —> y, and we can set everything in the 
integrands for K± constant, except — y. Hence 


K{Q,y) 

and 


\la-2 

Zmax. 


Vv max 


dz 

Vx-y’ 


K{7vl2,y) 


V2 


lla-2 

^max 


Vv max 


dz 

Vx-y' 


Thus 
r* 4 

Cy > y ^ ?/max. 


This is the minimum value of Cy. 


(56) 


(57) 


(58) 


Coefficients of Fricke expansion 


We derive an expression for Q„{9) (equation 44). First we 
define 


y = k(y+ + y-) 

^ = |(y+-?/-), 

and then set 


(59) 


y = y — Asin l3 


(60) 
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in equation (42). We write 


n p 

= / z"^'^~"{y+ + y-r~’^iy+-y-rdz. (6i) 

m=0 

Since g{y) is even in y we restrict our attention to even n. 
In this case the result is 


7r/2 


= J 2 I-" {^) f sin™/3 d/3, for even m; 


otherwise. 


Using standard integrals: 

r(i^) 


(”) 


2" r(i + f) 


(62) 


(63) 


for even m. In principle the integral in equation (61) is also 
analytic, since (n, m) are integers. 


This paper has been produced using the Royal Astronomical 
Society/Blackwell Science TJ/jX macros. 
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